Example: Blowfly Model
Likelihood-free inference for the blow-fly model was introduced by Simon N. Wood. We model here the discrete time stochastic dynamics of the size $N$ of an adult blowfly population as given in section 1.2.3 of the supplementary information.
where $e_t$ and $\epsilon_t$ are independent Gamma random deviates with mean 1 and variance $\sigma_p^2$ and $\sigma_d^2$, respectively.
using Distributions, StatsBase, LikelihoodfreeInference
Base.@kwdef struct BlowFlyModel
burnin::Int = 50
T::Int = 1000
end
function (m::BlowFlyModel)(P, N₀, σd, σp, τ, δ)
p1 = Gamma(1/σp^2, σp^2)
p2 = Gamma(1/σd^2, σd^2)
T = m.T + m.burnin + τ
N = fill(180., T)
for t in τ+1:T-1
N[t+1] = P * N[t-τ] * exp(-N[t-τ]/N₀)*rand(p1) + N[t]*exp(-δ*rand(p2))
end
N[end-m.T+1:end]
endLet us plot four realizations from this model with the same parameters.
using StatsPlots
plotly()
m = BlowFlyModel()
plot([plot(m(29, 260, .6, .3, 7, .2),
xlabel = "t", ylabel = "N", legend = false) for _ in 1:4]...,
layout = (2, 2))
To compare different realizations we will use histogram summary statistics. In the literature one finds also other summary statistics for this data.
summary_statistics(N) = fit(Histogram, N, 140:16:16140).weightssummary_statistics (generic function with 1 method)We will use a normal prior on log-transformed parameters.
function parameter(logparams)
lP, lN₀, lσd, lσp, lτ, lδ = logparams
(P = round(exp(2 + 2lP)),
N₀ = round(exp(4 + .5lN₀)),
σd = exp(-.5 + lσd),
σp = exp(-.5 + lσp),
τ = round(Int, max(1, min(500, exp(2 + lτ)))),
δ = exp(-1 + .4lδ))
end
(m::BlowFlyModel)(logparams) = m(parameter(logparams)...)
target(m::BlowFlyModel) = [(log(29) - 2)/2,
(log(260) - 4)*2,
log(.6) + .5,
log(.3) + .5,
log(7) - 2,
(log(.2) + 1)/.4]
lower(m::BlowFlyModel) = fill(-5., 6)
upper(m::BlowFlyModel) = fill(5., 6)
prior = TruncatedMultivariateNormal(zeros(6), ones(6),
lower = lower(m), upper = upper(m))TruncatedMultivariateNormal{Distributions.MvNormal{Float64,PDMats.PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},Float64}(
mvnormal: DiagNormal(
dim: 6
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Σ: [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]
)
lower: [-5.0, -5.0, -5.0, -5.0, -5.0, -5.0]
upper: [5.0, 5.0, 5.0, 5.0, 5.0, 5.0]
)
Let us now generate some target data.
model = BlowFlyModel()
x0 = target(model)
data = summary_statistics(model(x0))1000-element Array{Int64,1}:
0
0
0
0
0
1
0
2
0
5
⋮
0
0
0
0
0
0
0
0
0Adaptive SMC
smc = AdaptiveSMC(prior = prior)
result = run!(smc, x -> summary_statistics(model(x)), data,
maxfevals = 2*10^5, verbose = false)
using PrettyTables
pretty_table([[keys(parameter(zeros(6)))...] quantile(smc, .05) median(smc) mean(smc) x0 quantile(smc, .95)],
["names", "5%", "median", "mean", "actual", "95%"],
formatter = ft_printf("%10.3f"))┌───────┬────────────┬────────────┬────────────┬────────────┬────────────┐
│ names │ 5% │ median │ mean │ actual │ 95% │
├───────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│ P │ -2.192 │ 0.472 │ 0.046 │ 0.684 │ 2.078 │
│ N₀ │ -1.812 │ -0.109 │ -0.028 │ 3.121 │ 1.822 │
│ σd │ -1.740 │ -0.008 │ -0.065 │ -0.011 │ 1.413 │
│ σp │ -1.680 │ -0.054 │ 0.004 │ -0.704 │ 1.947 │
│ τ │ -1.479 │ 0.291 │ 0.206 │ -0.054 │ 1.724 │
│ δ │ -1.490 │ 0.124 │ 0.104 │ -1.524 │ 1.677 │
└───────┴────────────┴────────────┴────────────┴────────────┴────────────┘histogram(smc)
corrplot(smc)
KernelABC
k = KernelABC(prior = prior, delta = 1e-1, K = 10^3, kernel = Kernel())
result = run!(k, x -> summary_statistics(model(x)), data)
pretty_table([[keys(parameter(zeros(6)))...] quantile(k, .05) median(k) mean(k) x0 quantile(k, .95)],
["names", "5%", "median", "mean", "actual", "95%"],
formatter = ft_printf("%10.3f"))┌───────┬────────────┬────────────┬────────────┬────────────┬────────────┐
│ names │ 5% │ median │ mean │ actual │ 95% │
├───────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│ P │ -1.756 │ 0.541 │ 0.217 │ 0.684 │ 2.059 │
│ N₀ │ -1.563 │ 0.035 │ 0.072 │ 3.121 │ 1.784 │
│ σd │ -1.444 │ 0.052 │ 0.055 │ -0.011 │ 1.732 │
│ σp │ -1.714 │ -0.057 │ -0.030 │ -0.704 │ 1.794 │
│ τ │ -1.535 │ 0.155 │ 0.102 │ -0.054 │ 1.728 │
│ δ │ -1.456 │ 0.050 │ 0.092 │ -1.524 │ 1.770 │
└───────┴────────────┴────────────┴────────────┴────────────┴────────────┘histogram(k)
Kernel Recursive ABC (with callback)
k = KernelRecursiveABC(prior = prior,
K = 100,
delta = 1e-3,
kernel = Kernel(bandwidth = Bandwidth(heuristic = MedianHeuristic(2^3))),
kernelx = Kernel());We will use a callback here to show how the estimated parameters evolves.
using LinearAlgebra
res_krabc = run!(k, x -> summary_statistics(model(x)), data,
maxfevals = 1300,
verbose = true,
callback = () -> @show norm(k.theta - x0)/norm(x0))(x = [1.0291419932049777, 0.38359093735586475, 0.13204852686825827, 0.02008541956459824, -0.019452629431486828, 0.2492449358797572],)This page was generated using Literate.jl.